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Abstract 

An edge-based finite element formulation with vector absorbing 
boundary conditions is presented for scattering by composite struc- 
tures having boundaries satisfying impedance and/or transition con- 
ditions, Remarkably accurate results are obtained by placing the mesh 
a small fraction of a wavelength away from the scatterer. 
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1 Introduction 


Of generic interest in electromagnetic scattering is the modelling of composite 
configurations comprised of metallic and non-metallic sections. In the case 
of man-made structures, abrupt material discontinuities and metallic cor- 
ners are also encountered along with resistive sheets and thin ferrite coatings 
intended for controlling the scatterer’s radar cross-section (RCS). Differen- 
tial equation methods, especially the finite element method (FEM), with its 
capability of handling arbitrary geometries and its versatility in modelling 
inhomogeneities and material discontinuities has been a viable solution ap- 
proach for bounded domain problems. However, for unbounded problems 
as is the case with electromagnetic scattering, the solution is more involved 
since the finite element mesh needs to be truncated artificially at some dis- 
tance from the object with a suitable boundary condition. These boundary 
conditions can be either global or local. Global boundary conditions are ex- 
act but lead to fully populated submatrices thus spoiling the sparse, banded 
structure of the finite element system. In contrast, local conditions such as 
the absorbing boundary conditions( ABCs), are approximate but have the im- 
portant advantage of retaining the sparsity of the matrix system. ABCs are 
essentially differential equations enforced at the mesh truncation boundary 
and are chosen to suppress all non-physical reflections from that boundary, 
thus ensuring the outgoing nature of the waves. 

A variety of ABCs have been derived and widely employed in FEM so- 
lutions of open region two-dimensional scattering problems. However, the 
method’s implementation and performance for scattering by three dimen- 
sional geometries using edge-based finite elements has not received simi- 
lar attention. The only three-dimensional implementations of the FEM for 
scattering has been a hybrid solution combined with the boundary element 
method (BEM)[1,2] and a node-based formulation combined with ABCs[3]. 
The boundary element method, though exact, is equivalent to employing a 
global boundary condition for terminating the mesh and consequently leads 
to a full matrix, restricting the method’s utility to small geometries. For 
large-scale three-dimensional applications, it is necessary to employ an ABC 
for terminating the mesh to retain the 0(N ) storage requirement, character- 
istic of the finite element method. However, the use of traditional node-based 
elements for solving electromagnetic problems complicates the enforcement of 
the required boundary conditions at material discontinuities and impedance 
boundaries or resistive transitions. Moreover, the inherent singularities of 
the electromagnetic fields at metallic corners limits their use altogether. 

To avoid these difficulties, we consider an implementation of the FEM 
using vector basis functions whose degrees of freedom are associated with 
the fields along the six edges of a tetrahedron. In contrast to the traditional 
node-based elements, edge elements can treat geometries with sharp edges 
and are divergenceless. Moreover, they exhibit tangential continuity and 
normal discontinuity across inter-element boundaries and material disconti- 
nuities, thus making them ideal for representing electromagnetic fields. They 
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do, however, lead to a higher number of unknowns but this is balanced by 
the greater sparsity of the resulting finite element matrix. Thus the compu- 
tation time required to solve such a system iteratively, with a given degree 
of accuracy, is less than the traditional node-based approach[4]. Our imple- 
mentation is further coupled with a mesh termination scheme based on the 
vector ABCs derived in [5]. In contrast to the implementation proposed in 
[6], the one presented here preserves the symmetry of the finite element sys- 
tem, thus nearly halving the storage requirement and making it ideally suited 
for solution via a conjugate gradient type of algorithm. Further, the edge- 
based implementation discussed in [6] requires that the absorbing boundary 
be placed nearly a wavelength away from the scatterer, whereas in our imple- 
mentation remarkably accurate results are obtained with the ABCs enforced 
a small fraction of a wavelength from the scattering body. 


2 Formulation 

Let us consider the problem of scattering by an inhomogeneous target asso- 
ciated with possible material discontinuities. To solve for the scattered fields 
via the FEM, it is necessary to enclose the scatterer- embedded inside the 
volume V- by an artificial surface S 0 on which the ABC is enforced (see figure 
1). The ABCs to be considered in this paper are the Sommerfeld radiation 
condition given by 


n x V x E s = — jl 0 nxnxE’ (1) 

and the second-order ABC which can be written as 

n x V x E s = aE; + ^Vx [n(V x E*)„] + 0V,(V • E*) (2) 

where a = jk,/3 = l/(2jk + 2/r), E s represents the scattered electric field, 
n is the unit normal to the surface S 0 and the subscripts t and n denote 
the transverse and normal component to S a , respectively. When these ABCs 
are employed on the artificial boundary S 0 , they annihilate all field terms of 
0(r - ( 2m+1 )) and smaller, where m denotes the order of the ABC. The ABCs 
outlined above were derived for spherical surfaces but in this work we have 
extended their application to S 0 which include flat sections. This permits 
the construction of termination boundaries conformal to the scatterer, thus 
reducing the size of the the computational domain. 

The vector ABCs (1) and (2) can be combined and more conveniently 
written as 


n x V x E s = P( E s ) (3) 

for the scattered field formulation in which E s is the working variable and 

lixVxE = P( E) + U ,nc (4) 
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for the total field formulation where the unknown is the total electric field • 

In (4), 

u ,nc = n X V x E inc - P (E inc ) (5) 


where E = E s + E ,nc is the total field and E mc is the incident electric field. 
Considering (4) to be the boundary condition employed at S 0 , we can express 
the functional for the total electric field as 


F(E) 


[ [—(V x E) • (V x E) - 

JV \_f^r 

+ I [E • P( E) + 2 E • U ,nc ] 

J Sq 


k 2 a e r EE 

dS 


dV 


( 6 ) 


where t T and \i T are the relative permittivity and permeability, respectively. 

The above functional can be generalized to account for the presence of 
impedance and resistive sheets or other dicontinuous boundaries. In the case 
of a resistive card, the transition condition[7] 


n x (n x E) = -Rh x (H + - H") (7) 

must be enforced, where H* denotes the total magnetic field above and below 
the sheet, R is the resistivity in Ohms per square and n is the unit normal 
to the sheet pointing in the upward direction (+ side). For an impenetrable 
impedance surface, the appropriate boundary condition on that surface is 


n x (n x E) = —rjh x H 


( 8 ) 


where n is the unit normal to the surface and 77 is the surface impedance. Tak- 
ing into consideration these boundary /transition conditions, the functional 
for the total electric field can be more explicitly written as 


F( E) = / 

JV 


v in 


-(V x E) • (V x E) - k 2 0 e T E • E 


dV 


+jk 0 Z o [ -^7 (n x E) • (n x E) dS 
Js k A 

+ [ [E • P(E) + 2 E • U inc ] dS 

J Sq 


( 9 ) 


where K is the surface resistivity(i?) when integrating over a resistive card 
and equals the surface impedance^) for an impedance sheet. 

The formulation presented above is in terms of the total field but we can 
easily revert to a scattered field formulation by setting E 4 = E — E mc and 
noting that the scattered field satisfies the wave equation inside the domain 
of interest. The functional F(E S ) is given by 


f < E *> = L 


— (V x E s ) • (V x E s ) - k 2 0 e r E a • E 5 

hr 


dV 
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+jk o z o [ ~7 (n x E 3 ) • (n x E 3 ) dS 
Js k A 

+ / E 3 • P(E 3 )d5 

Js 0 

+2 jk 0 Z 0 [ — E 3 • (n x H ,nc ) dS 

JS d Hr V ' 

2 l 


+2 


(V x E 3 ) • (V x E ,nc ) - k 2 0 e r E 3 • E 


dV 


+2 jk 0 Z 0 f l(n x E 3 ) • (n x E inc )dS 
J Sf- F 

+/ (E'“') 


( 10 ) 


where Vd is the volume occupied by the dielectric (portion of V where e r or 
Hr are not unity), Sd encompasses all dielectric interface surfaces and 

Js E s • P(E s )dS = ^ [a (E 3 ) 2 + /?(V x E 3 )^ - /9(V • E 3 ) 2 ] dS 

when the second order ABC is employed. The function / (E ,nc ) is solely 
in terms of the incident electric field and vanishes when we take the first 
variation of F(E 3 ). We remark that the scattered field formulation was 
implemented in our code; however, we expect that the total field formulation 
would yield comparable results. 


3 Finite element discretization 

To discretize the functional given in (10), the volume V is subdivided into a 
number of small tetrahedra, each occupying the volume V e (e = 1, 2, ■ • • , M), 
where M denotes the total number of tetrahedral elements. Within each 
element, the scattered electric field is expressed as 

771 

E e = = {W e } r {F e } = {£ e } T {W e } (11) 

3 = 1 

where are the edge-based vector basis functions[4], E? denote the expan- 
sion coefficients of the basis and represent the field components tangential 
to the jth edge of the eth element, m is the number of edges making up 
the element and the superscript stands for the element number. The basis 
functions used in our implementation have zero divergence and constant curl. 

The system of equations to be solved for Ej is obtained by a Rayleigh- 
Ritz procedure which amounts to differentiating F with respect to each edge 
field and then setting it to zero. On substituting (11) into (10), taking the 
first variation in F and assembling all M elements, we obtain the following 
augmented system of equations 

( 1 M A/ * A/ p 

{aF } = + D B ’K £ '} + Efc") = o (i2) 
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In this, M s denotes the number of triangular surface elements on Sk and S a 
and M p is equal to the sum of the surface elements on S'*, Sd and the volume 
elements in Vj. The elements of the matrices [A e ], [i? s ] and {C p } are given 
by 


A‘j = [ [— (V x W-) • (V x Wj) — l’ 0 2 £ r W' • 

J JV e [fi r 

SI = iKZ, U jr(A x WJ) • (n x W‘)dS 

+ JJa Wf, ■ W-, + /S(V x W*)„ • (V x W*)„ - 0(V ■ WJ)(V • w;,)]<JS 

Cf = 2 jhzAl — Wf- (n x H ,W U? + / 4(n X w 0 • (ft x E"‘)iS 
[ J S P d Hr V 1 JS p K J 


+2 / —(V x Wf) • (V x E inc ) - k*e T W p ■ E ,nc 

Jv S Mr J 


dV 


where VJ is the volume of the pth tetrahedron inside the dielectric, S 3 and 
S p represent the surface area of the sth and pth triangular surface element 
and the subscripts t and n denote the tangential and normal components of 
a vector, respectively. The boundary condition n x E' = — n x E ,nc must 
be imposed a priori on metallic boundaries; however, no special treatment is 
required at material discontinuities. Only the identification of the edges on 
material discontinuities or inhomogeneities is required to kick in the contri- 
bution from the surface integrals in (10). 

The biconjugate gradient algorithm was used to solve the sparse, sym- 
metric system of equations. The residual norm was usually set to less than 
0.1% of the solution norm as a criterion for convergence since lower tolerances 
did not appear to offer significant improvement on the far-field values. The 
data structure was constructed such that only the non-zero elements of the 
upper triangular part of the symmetric, sparse matrix were stored in a iV a x k 
complex array. In our case, N a was typically 1.1 x N u , where N u denotes the 
number of unknowns and k was equal to 12. The corresponding addresses 
were stored in a separate N a x k integer array. The storage required in this 
scheme was about 25 N u and the number of distinct non-zero elements was 
typically 9N U . 


4 Results 

A computer program was written for implementing the proposed FE-ABC 
formulation. This implementation was validated by computing the scattering 
for several configurations including metallic and dielectric bodies a s well as 
structures satisfying resistive and impedance boundary conditions. 

Figure 2 compares the measured[8] bistatic cross-section (0 ,nc = 180°, 
<f> tnc = 90°) of a metallic cube having an edge length of 0.755A with the 
corresponding pattern computed by the three-dimensional FE-ABC code. 
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The second-order vector ABC was employed on a spherical mesh truncation 
boundary which was placed only 0.1 A from the edge of the cube. About 
33,000 unknowns were used for the discretization of the computational do- 
main and the [A] matrix contained a total of 264,000 distinct non-zero entries. 
The storage requirement of this matrix was consequently much smaller than 
that of the 1400 unknown moment method system (assuming the same sam- 
pling rate as the FEM of 14 points/A) which had 2 million non-zero entries. 

In figure 3, we plot the normal incidence backscatter RCS of a perfectly 
conducting cube as a function of its edge length. The meshes constructed 
for this experiment were terminated on conformal boundaries, i.e, on an- 
other cube placed a small distance (more than 0.15A) from the scatterer. As 
seen, the agreement with measured data[9] is remarkably good over a 50dB 
dynamic range. 

Figure 4 presents backscatter data for a cylinder of radius 0.3A and height 
0.6A. The data from the three-dimensional finite element code again compare 
well with that obtained from a moment method-body of revolution code. The 
mesh was terminated on a spherical boundary at a distance of 0.3A from the 
edge of the scatterer and the system consisted of nearly 33,000 unknowns. 
Convergence was achieved within about 350 iterations when the Sommerfeld 
radiation condition was imposed on the spherical mesh termination boundary. 
Each iteration took approximately 0.3 seconds on a Cray YMP and on the 
average it was found that for N > 25, 000, the number of required iterations 
were approximately N/100. The agreement was quite good even on enclosing 
the metallic cylinder in figure 3 with a rectangular outer boundary placed 
0.3A from the edge of the scatterer. 

The results presented till now have been for perfectly conducting ge- 
ometries. However, the real advantage of the FEM over integral equation 
techniques is the ease with which the former can handle material inhomo- 
geneities and transition conditions. With this in mind, the remaining figures 
show backscatter and bistatic patterns for scatterers comprised of resistive 
cards, dielectric material and combinations of these. One of the test cases 
was a prolate spheroid shown in figure 5 filled with lossy dielectric having 
a permittivity of 4 — j 1, k 0 a = tt/2 and a/b = 2, where a and b are the 
major and minor axes of the spheroid, respectively. The bistatic pattern 
(0 tnc = 180°; <p' nc = 90°) obtained from the FE-ABC solution agree reason- 
ably well with those obtained via the hybrid finite element-boundary integral 
method presented in [2]. However, the corresponding convergence rate for 
non-metallic bodies and resistive/impedance sheets was found to be slower 
than that observed for metallic scatterers. A diagonal preconditioner was, 
therefore, used to accelerate the convergence of the biconjugate gradient al- 
gorithm with encouraging results. 

For our last example, we compute the scattering from an inhomogeneous 
geometry with embedded resistive cards. Particularly, the scatterer shown in 
figure 6 consists of an air-filled resistive card block (0.5Ax0.5Ax0.25A) joined 
to a metallic block (0.5A x 0.5A x 0.25A). In figure 7, we compare a principal 
plane backscatter pattern obtained from our 3D FE-ABC implementation 
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with data computed using a traditional moment method code[10]. For the 
FE-ABC solution, the scatterer was enclosed within a cubical outer boundary 
placed only 0.3A away from the scatterer. This resulted in a 30,000 unknown 
system which converged to the solution in about 400 iterations when using 
the Sommerfeld radiation condition and in 1600 iterations when the second 
order ABC was used. For this geometry, the second order ABC did not 
provide a significant improvement in accuracy (only about O.ldB) over the 
first order condition. The same case was run with a higher discretization 
resulting in a system of 50,000 unknowns; however, there was no significant 
difference in the far-held values with the earlier case. The geometry for the 
backscatter pattern shown in figure 8 is the same as in figure 6 with the 
air- filled section now occupied by a lossy dielectric having e T = 2 — j2. The 
backscatter echo-area pattern for the <f><i> polarization as computed by our 
FE-ABC code is again seen to be in good agreement with corresponding 
moment method data[10]. 


5 Conclusions 

In this paper, we have shown that the finite element technique with vector 
basis functions, when coupled with ABCs for mesh termination and the bi- 
conjugate gradient algorithm for the solution of the resulting system, is a 
viable procedure for computing the scattering by three-dimensional targets. 
We have found that these ABCs can be enforced only a small fraction of 
a wavelength from the scatterer’s surface. This is probably due to the fast 
(1/r) decay of the scattered fields. As a result, in addition to the sparsity 
of the matrix, the total number of unknowns is kept under control. Fur- 
ther, due to the use of edge elements, the program can easily handle sharp 
conducting edges and tips, inhomogeneous dielectric and/or magnetic ma- 
terials, resistive sheets and impedance surfaces. These, in conjunction with 
the well-known advantages of the finite element method, results in low O(N) 
storage requirement, making the computation of large body scattering possi- 
ble. These capabilities along with the ease in modelling arbitrary geometries, 
makes this formulation, to the best of our knowledge, one of the first suitable 
for solving practical three-dimensional scattering problems. 

6 Appendix 

6.1 Anisotropic case 

In order to deal with anisotropic scatterers, the functional outlined in (9) 
undergoes a slight modification since the material properties of the scatterer 
(permeability and permittivity) are now second rank tensors rather than 
scalars. Equation (9) can therefore be written as 

n E) = / y [(VxE).{p3-'.(VxE)}-ljE.{f.B)]iV 
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( 13 ) 


+jk 0 Z 0 f (n x E) • (n x E) dS 
J Sk A 

+ J s [e • 5(E) + 2E • U’ nc ] dS 

where 

fixx fixy fixz 

fiyx fiyy fiyz ( 14 ) 

fizx fizy fizz 

and 

€xx ^xy ^xz 

C r = €yx Cyy Cy Z 

t zx tzy € zz 

The changes in the scattered field functional (10) and in the elements of the 
matrices [A e ], [5 s ] and { C p } can also be easily incorporated by replacing 

— X-Y with X-t^T'-Y 

/*r 

and 

e r X • Y with X-e^-Y 
where X and Y are arbitrary vectors. 
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Figure captions 


Figure 1: Illustration of a scatterer enclosed by an artificial surface, S a , on which 
the absorbing boundary condition is imposed. 

Figure 2: Bistatic pattern of a metallic cube having an edge length of 0.755A. 
The dotted circles and squares are measured data. 

Figure 3: Backscatter RCS at normal incidence for a metallic cube of side a. The 
dots indicate measured data. 

Figure 4: Backscatter pattern of a perfectly conducting cylinder of radius 0.3A 
and height 0.6A. The axis of the cylinder coincides with the z-axis of 
the cartesian coordinate system. The dotted circles and squares are 
data obtained from a moment method code. 

Figure 5: Bistatic pattern of a lossy prolate spheroid (e r = 4 — j 1; k 0 a = 7t/2; 

a/b = 2), where a and b are the major and minor axis of the spheroid, 
respectively. The dotted circles are data computed from a hybrid 
FEM/MoM solution. 

Figure 6: Geometry of a cube (a=b=0.5A) consisting of a metallic section and 
an air-filled section, where the latter is bounded by a resistive surface 
having R = Z D . The outer cube drawn with dashed lines represents the 
surface on which the ABC was enforced. 

Figure 7: RCS patterns in the x-z plane for the geometry in figure 6 with H z nc = 
0(solid curve is the FE-ABC pattern; circles are MOM data) and E' z nc = 
0(dashed curve is the FE-ABC pattern; black dots are MOM data). 

Figure 8: RCS pattern in the x-z plane for the geometry in figure 6 with the 
air-filled section replaced with dielectric having e T = 2 — j 2. The solid 
curve is the FE-ABC pattern and the black dots are MoM data for the 
E' z nc = 0 polarization. 
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Figure 1. Illustration of a scatterer enclosed by an artificial 
surface, So , on which the absorbing boundary 
condition is imposed. 
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Figure 2: Bistatic pattern of a metallic cube having an edge length of 0.755A. 
The dotted circles and squares are measured data. 
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Figure 4: Backscatter pattern of a perfectly conducting cylinder of radius 0.3A 
and height 0.6A. The axis of the cylinder coincides with the z-axis of 
the cartesian coordinate system. The dotted circles and squares are 
data obtained from a moment method code. 
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Bistatic pattern of a lossy prolate spheroid (e r = 4 “ J }' k ° a ~ *( ' 
alb = 2), where a and b are the major and minor axis of the spheroid, 
respectively. The dotted circles are data computed from a hybrid 
FEM/MoM solution. 
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Figure 6: Geometry of a cube (a=b=0.5A) consisting of a metallic section and 
an air-filled section, where the latter is bounded by a resistive surface 
having R = Z a . The outer cube drawn with dashed lines represents the 
surface on which the ABC was enforced. 



Figure 7: RCS patterns in the x-z plane for the geometry in figure 6 with H' z ne 
0(solid curve is the F E-ABC pattern; circles are MOM data) and E\ nc 
Ofdashed curve is the FE-ABC pattern; black dots are MOM data). 




Figure 8: RCS pattern in the x-z plane for the geometry in figure 6 with the 
air-filled section replaced with dielectric having e r = 2 — j2. The solid 
curve is the FE-ABC pattern and the black dots are MoM data for the 
E\ nc = 0 polarization. 



